In [2]:
from __future__ import division
from scipy.stats import norm
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import prettyplotlib as ppl
import brewer2mpl
mpl.rcParams['lines.linewidth'] = 2
mpl.rcParams['lines.color'] = 'r'
mpl.rcParams['axes.titlesize'] = 32
mpl.rcParams['axes.labelsize'] = 24
mpl.rcParams['axes.labelsize'] = 24
mpl.rcParams['xtick.labelsize'] = 24
mpl.rcParams['ytick.labelsize'] = 24
colors = brewer2mpl.get_map('Set1', 'qualitative', 9).mpl_colors
def plot_power_curve(alpha, stddev, sample_sizes, alternative=">"):
if alternative == ">":
z_crit = norm.interval(1 - alpha * 2)[1]
x = np.arange(0, 14, 0.01)
elif alternative == "<":
x = - np.arange(0, 14, 0.01)
z_crit = norm.interval(1 - alpha * 2)[0]
plt.figure(figsize=(12, 8))
for i, n in enumerate(sample_sizes):
z = (-x / (stddev / np.sqrt(n))) + z_crit
z = -z if alternative == "<" else z
ppl.plot(x, 1 - norm.cdf(z), color=colors[i], linewidth=3, label="n = %d" % n)
plt.legend(loc=4, fontsize=24)
plt.xlabel("Difference ($\mu_{hyp} - \mu_{true}$)")
plt.ylabel("Power")
plt.title("Power curve for $\sigma = %d$ and $\\alpha = %.2f$" % (stddev, alpha))
plt.tight_layout()
plt.ylim([0, 1.1])
return plt
In [4]:
plt = plot_power_curve(0.05, 15, [20, 30, 50, 100, 1000])
In [19]:
x = np.arange(80, 200, 1)
plt.plot(x, norm.pdf(x, 160, 34/np.sqrt(10)), color=colors[0], label="hypothesized")
plt.plot(x, norm.pdf(x, 110, 34/np.sqrt(10)), color=colors[1], label="true")
plt.legend()
x = np.arange(140.9, 180)
plt.fill_between(x, norm.pdf(x, 110, 34/np.sqrt(10)), color='grey')
Out[19]: